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Melting in Monolayers : Hexatic and Fluid Phases. 
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There are strong evidences that the melting in two dimensions depends crucially on the form 
and range of the interaction potentials between particles. We study with Monte Carlo simulations 
the phase diagram and the melting of a monolayer of point-particles interacting with repulsive 
Inverse Power Law Interactions, V(r) = Q 2 (a/r) n where n can take any real positive value (n-OCP 
monolayer). As n is varied from to oo (Hard Disks), including Coulomb (n = 1) and Dipolar 
(n — 3), melting occurs with different mechanisms and the overall picture permits to understand 
the diversity of mechanisms found experimentally or in computer simulations for 2D melting. The 
empirical transition curves for n < 3 and the excellent qualitative and semi-quantitative agreements 
with the KTHNY theory found for the melting of n-OCP monolayers with n < 3 are the main 
results of the present work. 



PACS numbers: 64.70.D-, 61.20.Ja, 64.60.De 



I. INTRODUCTION 



The crystalline phases in two dimensional systems are 
characterised by long-range orientational order and, be- 
cause of the importance of fluctuations 0, Hj], a quasi- 
long range positional order rather than a long ranged 
order. Several experimental systems, as electrons at 
surface of liquid Helium KMi or in MOSFETs transi- 
tors (Wigner crystals) [!, uHlO] , dusty plasma [lil ] , col- 
loids confined in slabs or at interfaces |12| - H7| , Abrikosov 
vortex in 2D supraconductors_ 191. 20| , rubidium atoms 



trapped in 2D optical traps [21, 22 j, atoms adsorbed at 
liquid or solid surfaces [23j], may be considered as 2D 
systems. Most of these systems exhibit phase transi- 
tions between their crystalline and disordered (hexatic 
or isotropic) fluid-like phases. To describe 2D meltings, 
a variety of mechanisms has been proposed [241 ] . am ong 
them, there are the grain-boundary induced melting [25] 
where the Solid-Fluid transition is described as a first 
order phase transition between the solid and the fluid 
phases, and the KTHNY theory (after Kosterlitz, Thou- 
less, Halperin, Nelson and Young) [26-33], with an inter- 
mediate hexatic phase and two continuous phase transi- 
tions (Solid/Hexatic and Hexatic/Fluid). The melting in 
two dimensions may also be described in the vicinity of 
the isostructural critical points (3J-|36[ . The isostructural 
transition describes a Solid-Solid transition in which the 
two Solid have the same symmetry but a different lattice 
spacing. The hexatic phases in systems with isostructural 
transition occur when the triple point (Fluid-Solid-Solid) 
is close enough to the isostructural critical point (Solid- 
Solid) [MSI- 

From experimental and numerical simulations results, it 
is clear that the underlying mechanisms for 2D meltings 
depend on the geometry of the systems [13, [3 and inter- 



actions between constituants ; in particular, the melting 
depends strongly on the range and form of the interac- 
tion pa, m m m mm- 

In this paper, we report a preliminary study of the melt- 
ing in monolayers of point particles that interact with 
repulsive Inverse Power Law (IPL) potentials. The re- 
pulsive Inverse Power Law interactions are 



(1) 
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where Q 2 is the energy scale, a the length scale and r 
the distance between the particles. In this class of po- 
tential, the power n can be considered as a parameter 
and it shall take any real positive value in the follow- 
ing [49-56]. Although n is considered in this paper as 
a continuous parameter, some integer values of n corre- 
spond to physical interactions in experimental and real 
systems, as for instance : Coulomb interactions (n = 1) 
[lO, El, S3> dipolar interactions (n = 3) [ll |5S|. 
the short range repulsion in the Lennard- Jones poten- 
tial (n = 12) [l8|, |43j-|46j and also the hard disk interac- 
tion in the limit n — ► oo (then a is the diameter of the 
disks) [38H41 Hfl ■ For all these interactions the melt- 
ing in two dimensions has been studied experimentally 
and numerically ; different mechanisms have been found 
including the grain boundary induced melting [24l l25"j , 
the KTHNY theory [H [US, S El and also, on hard 
disks systems, a first order transition from the isotropic 
fluid to the hexatic phase followed by a continuous transi- 
tion from the hexatic phase to the 2D-hexagonal crystal 
38]. Since the interaction between the particles in all 
these systems belongs to the class of repulsive Inverse 
Power Law potentials, by considering the power n as a 
continuous parameter we are able to define a model on 
which the crossover between all the mechanisms for the 
two dimensional melting can be studied. 
It is worthwhile to outline also that repulsive Inverse 
Power Law Interactions are used as repulsive reference 
potentials in dense liquids to study the thermodynamic 
and dynamical properties [54 56]. In these models, IPL 
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potentials are effective interactions and the exponent n 
is computed from the scaling properties of the diffusion 
coefficients (EHJ or by using the virial theorem [HH ; the 
values of n computed likewise are in general greater than 
fO. Repulsive Inverse Power Law interactions are also 
used in studies on glass transitions [H, l59U6ll ] . 
The purpose of this paper is to define the n-OCP mono- 
layer model and to give a first overview of its properties 
and phase diagram, most of the results reported in the 
present work are obtained for long ranged and quasi long 
ranged interactions (n < 2 — 3). 

The paper is organised as follows. In section II, we de- 
scribe the numerical methods and the model system used 
in this study and we give computational details on the 
Monte Carlo simulations. Based on the results obtained 
by simulations, the empirical transition curves as func- 
tion of n are given in section III. In section IV, we re- 
port the results of Monte Carlo simulations done in the 
present work. In this section, a particular attention is 
paid to the Kosterlitz-Thouless essential singularity [3l[ 
near the Fluid/Hexatic transition, since this singularity 
is a clear signature of the KTHNY theory 26] . A discus- 
sion of the results reported in the paper is done in section 
V. 

For completeness, in an Appendix we give the Ewald 
sums for repulsive Inverse Power Law interactions in 
monolayers, a complete analytical derivation can be 
found in refs. 14ft 5(1 



II. 



METHODS AND COMPUTATIONAL 
DETAILS 



To explore the phase diagram of a system of point par- 
ticles with pairwise repulsive inverse power law interac- 
tions we use Monte Carlo simulations. The computations 
are done in the canonical ensemble with a variable shape 
of the simulation box, but at a constant surface area A 
[67l l68j ] . The Ewald method for inverse power law inter- 
actions derived in refs.[4{||5(| is used in all computations. 
For long and intermediate ranged potentials (n < 7) the 
use of Ewald methods is necessary to avoid bias in compu- 
tations. If the interactions are very long ranged (n < 2), 
a neutralizing background is used to obtain convergent 
sums as in the OCP model [69j ; in the following, the sys- 
tems of point particles with Inverse Power Law potential 
interactions Eq. (TTJ with a neutralizing background when 
convergence is needed [lij [5(| [z3l are called n-OCP sys- 
tems. 

For short ranged interactions (repulsive soft spheres - 
n > 10) truncations of the interactions can be done 
safely, however, in this study we use Ewald methods also 
for these values of n to avoid any type of discontinuity 
that could be introduced by changing the method used 
to compute energies. In the appendix, one would find 
the energy of the monolayer computed with the Ewald 
method. 

All thermodynamic points are obtained on systems with 



N = 4096 particles, equilibration is done by making 
2 x 10 5 MC-cycles and the number of MC-cycles to com- 
pute averages are between 2 x I0 5 and 2.3 x I0 6 . With 
the code written for this study, the cpu-time for one MC- 
cycle is about one or two seconds ; however, despite this 
quite good efficiency, the computation time needed to 
obtain one thermodynamic point is about one or two 
months. Improvements of the source code with mesh 
methods (PME, P3M, etc.) [49] can not be done easily 
for non integer value of the power n, especially for the 
grid points fractional charge assignment step [49j . 
The reduced units are defined as follows : density p* — 
pa 2 ; coupling T = Q 2 /ksT = Q* 2 , all computations are 
done with Q = 14. For this system, the only relevant 
parameter is the coupling constant T = (npa 2 ) n ^ 2 T. For 
notational convenience, the asterisks will be dropped in 
the remaining of the paper. 

Bond orientational order parameters (BOOP), their his- 
tograms and their correlations are computed by using 
Voronoi constructions. It is worthwhile to note that 
Voronoi (or Delaunay) constructions to compute order 
parameters and a variable shape of the simulation box 
are necessary to study fluid-solid phase transitions in two 
dimensions to avoid systematic bias 62]. Voronoi con- 
struction permits to identify exactly the neighbors of all 
particles and variation of the shape of the simulation box 
is necessary to avoid irrelevant hysteresis loops induced 
by some inadequacies between the geometry of the box 
and the geometry of the cristal or hexatic phases, even 
for large and very large systems. 

The finite size analysis presented in section IV is achieved 
by dividing the system in smaller subsystems in a similar 
way as it is done in ref. t 6j|, the fourth-order cumulants 
are defined as 



U 6 (N) = I 



< <t>6 >N 



3 < 



(2) 



N 



with N the average number of particles in subsystems. 
A more detailed study using different system sizes and 
the multiple histogram reweighting method (65 . 66] is on- 
going. Some results obtained with the multiple histogram 
method are given in section IV for systems with 4096 and 
2025 particles and for n = 1.65 and p = 0.109. 
To obtain measurements of the correlation lengths and 
exponents of the algebraic decay of correlations in hexatic 
and fluid phases, the tails of correlation functions (i.e. 
for s > 10 a) are fitted as follows. The tails of bond 
orientational correlation functions in the Fluid phase are 
represented as 

ge(s) =< <f> 6 > 2 g +dpe-'MV 

+G (B) e - s ,UT) j Q 

and, in Hexatic and Solid phases, as 

,>;6(T) 

S/ \ Ai: 



(3) 



g 6 (s) =< 6 > 2 g +H 6 (-J J o( 2 ^— ) ( 4 ) 



with Jo(x) the Bessel function. Similar forms are used 
for center to center correlation functions, with parameter 
£(T), rj(T) and A, as 

g(s) = 1 + Go e- s/ « (T) J (2^) (5) 

in the Fluid and Hexatic phases, and as 

Jo (2tt~) (6) 



g(s) = l + G s {-) J„ 



in the Solid phase. 

The center to center correlation functions g pq {s) for pair 
of particles with p neighbors and q neighbors are com- 
puted from the Voronoi constructions and normalised to 
lass-> oo. In Solid and Hexatic phases, since the parti- 
cles with a number of neighbors different from 6 are quite 
rare, large statistical fluctuations in the g P q{s) functions 
are found for s > 10a and extremely strong correlations 
are found for s ~ 1/ 'y/Wp. 

The Kosterlitz-Thouless essential singularity near transi- 
tions, for a system with N particles, can be represented 

asMMm 



X( r ,N) = X (N)exp 



a(N) 



l7~7c(JV) | l 



(7) 



with 7c (ao = i/r c (iv). 

For MC computations done with n, p and Q fixed, one 
can replace 7 c (iV) by T C (N). All quantities exhibiting 
a KT-singularity (correlation lengths and quantities re- 
lated to them) are fitted with Eq.([7]) to obtain estimates 
of T C (N) or T C (N) [13J. For a given set of computations 
(at fixed n and p, or at fixed p and T), two temperatures 
or coupling constants are found ; we identify one with 
the Fluid/Hexatic transition (Tp/jj or Tp/jj) and the 
other with the Hexatic/Solid transition (T H / S or T H / S ). 
We locate the Fluid/Hexatic transitions from the KT- 
singularities found for : the finite-size analysis of < (j>e > 
and X 6 (N -> oo - cf. Inset FIGgJb)), £(T) and £ 6 (T) 
; and from £+(T) for the Hexatic/Solid transitions. The 
values found for T c (n) by doing these analysis are re- 
ported on TABLE U and the value of p c (n) deduced from 
r c (n) at T = 1 is also given on TABLE Q] and represented 
as symbols denoted MC on FIGOJa). As explained in the 
next section, the values of T c (n) permit also to construct 
empirical transition curves, for the two transitions, as 
function of the power n of the IPL potentials. 



III. THE EMPIRICAL TRANSITION CURVES 
AS FUNCTION OF n 

The main purposes of this paper are to obtain a first 
overview of the phase diagram of the monolayer n-OCP 
model as function of the parameter n and to locate ap- 
proximately the critical or multicritical points in the 
(n, p) plane or (n, T) plane. The coupling constant T is 
the only relevant parameter for the n-OCP system. Pre- 
vious experimental and numerical studies on Coulomb, 
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FIG. 1: (Color online) Phase diagram of monolayers (T = 1.0 
and Q = 14). (a) Representation of the empirical transition 
curves for the Fluid/Hexatic and Hexatic/Solid transitions 
(Eq.Q). Inset : enlargement of the region 0.6 < n < 3, 
locations of F and H are respectively np ~ 1.63 and pfo 2 ~ 
0.130 and n H 1.88 and p H ° 2 oi 0.143. Runs A are MC 
computations done at irpa 2 — 1.35 and Runs B at pa 2 = 0.3. 
(b) Bond orientational order parameters < <f)Q > for Runs 
A and B. (c) Binder cumulants for Runs A near the F/S 
transition. The red lines are obtained by straight line fit of 
the cumulant in the region of transition, they are crossing for 
n ~ 12.85 — 12.90. (d) Same as (c), but for Runs B ; red lines 
cross around n = 7.35. 



dipolar and hard disk systems show that the value of r c 
at transition depend on the interaction. Therefore, we 
may represent the phase diagram of the monolayer n- 
OCP model in the (n, F) plane with the help of curves 
r c (n) or equivalently, in the (n, p) plane as curves 



Pc{n) 



1_ (T c {n)k B T 



Q 2 



2/n 



1 



r c (n) 



2/n 



(8) 



At present, we are not able to describe a theory that 
will give an analytical derivation of r c or p c as function 
of n. However, the numerical and the experimental re- 
sults obtained is previous works d, [9l-[Tll [l3l . l38j and the 
Monte-Carlo computations done in the present work al- 
low us to give an empirical representation of the r c (n) 
and p c (n) curves. 

On TABLE [H we report some numerical evaluations of 
r c (n) and p c (n) for several systems including : Coulomb 
interactions (n = 1) 0, fol— Tlll| , dipolar interactions be- 
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Fluid/Hexatic 


Hexatic/Solid 




n 


pc{n) 


r c (n) 


n 


p c (n) 


T c (n) 


Transitions References 


Methods 


HD 






HD 












(n = oo) 


0.891-0.912 


oo 


(n = oo) 


0.917 


oo 


1 st jr\nd 

1 


[381 


Simulations (MC) 


0.66 ± 0.09 


1.35/tt 


216.4 ±9.5 


0.71 ±0.05 


1.35/tt 


217.9 ±6.0 


KTHNY 


This work 


Monte Carlo (A) 


0.71 


0.312 ±0.018 


194.7 ± 1.2 


0.71 


0.388 ±0.011 


210.2 ± 1.0 


KTHNY 


This work 


Monte Carlo 


0.72 ± 0.02 


0.3 


191.9 ±2.9 


0.86 ±0.02 


0.3 


191.1 ±5.0 


KTHNY 


This work 


Monte Carlo (B) 


1 


0.123 - 0.192 


137 ± 15 


1 


0.123 - 0.192 


137 ± 15 




[3-6] 


Experiments 


1 


0.119 


120 


1 


0.163 


140 




40] 


Quantum Monte Carlo 


1 


0.123 


122 


1 


0.127 


124 


KTHNY 


m 


Monte Carlo 


1.48 


0.132 ±0.016 


102.3 ±0.5 


1.48 


0.152 ±0.025 


113.5 ± 1.0 


KTHNY 


This work 


Monte Carlo 


1.65 


0.128 ±0.010 


92.9 ±0.3 


1.65 


0.138 ±0.028 


98.3 ± 1.0 


KTHNY 


This work 


Monte Carlo 


3+ 


0.141 


57.5 


3+ 


0.147 


61.5 


KTHNY 


[13] 


Experiments 


3* 


[1/2^3] 


69.1 ±0.5 


3* 


[l/2%/3] 


72.0 ± 1.0 


KTHNY 


[58] 


Molecular Dynamics 


3 


0.159 ± 0.008 


69.3 ±0.4 


3 


0.161 ±0.019 


70.7 ± 1.0 


KTHNY 


This work 


Monte Carlo 


Fluid/Solid 


n 


p c (n) 




r c (n) 






7.30-7.40 




0.3 


157.4 - 157.9 


l at order 


This work 


Monte Carlo (B) 




12.75-12.95 




1.35/vr 




1328 - 1368 


1 st order 


This work 


Monte Carlo (A) 



See text ; * The reduced units used in ref. [58fl are equivalent to Q 2 = 1 and the computations are done at the density 



indicated in the table, this density does not correspond to the critical density at T — 1 as in the other lines of the table. 

TABLE I: Estimations of densities and n-OCP coupling constants for the Fluid/Hexatic (F/H) and Hexatic/Solid (H/S) 
transitions for several values of n at T = 1.0. Data estimated in this work are extracted from the behaviour of the bond 
orientational correlation function ga(s) computed with Monte Carlo simulations and from the Kosterlitz-Thouless essential 
singularities. These data give the empirical transition curves as in Eq.© and for n < 3 with : r o r/ "> ~ -0.84 ; T\ ' ' ~ -1.42 
; F ( 2 F/H) ~ -0.16 and T { H/S) ~ -0.74 ; r[ H/s) ~ -1.45 ; F< H/S) ~ -0.24. 



tween superparamagnetic colloids (n = 3) |kl[58|, sim- 
ulations on Hard Disks systems (n — > oo) [38| | and MC 
computations done in the present work. These data are 
quite well represented by choosing T c (n) as 



In 



r c (n) 
f 



T ±rilnn±r 2 (lnn) 2 + T C 



(9) 



We assume the same dependence of r c on n for both 
Fluid/Hexatic (F/H) and Hexatic/Solid (H/S) transi- 
tions but with different value of the T^-parameters for 
each transition. For both transitions, Too may be ob- 
tained from the numerical computations done recently 
by Bernard and Krauth on hard disks systems [IK . The 
asymptotic behavior of Eq.@ 



2 i im ^r e (n)/f\ =ln 2 



(10) 



gives, with p, 



(F/H) 



JH/S) 



0.917 m 



t£ /H) ~ 0.515 and 



0.891 and p\ 

0.529. The transitions 
curves F/H and H/S obtained by using the empirical fits 
of Eq.© are represented on FIG [Ha). It is worthwhile 



to note that the transition curves as in Eq. ((SJ are equiv- 
alent to those found by Platzman and Fukuyama [H, 0] 
for the Coulomb OCP monolayers when the Linedeman 
criterium is applied in the classical region. 
On TABLE HI we interpret the differences found for n — 3 
between MC computations and the experiments done by 
P. Keim and co-workers as follows. In the very nice ex- 
periments done in refs. 13], the superparamagnetic col- 
loids are confined at an air/water interface and dipolcs 
are induced by an external magnetic field ; therefore a 
small fluctuation in the vertical position of the parti- 
cles would result in a small attractive contribution in 
the dipolar interactions, this would make the ordered 
phases (Hexatic and Solid) stable at smaller densities, 
smaller values of the coupling constant. This may ex- 
plain the small quantitative differences, but an excellent 
qualitative agreement, found between superparamagnetic 
colloids confined at an air/ water interface and pure nu- 
merical monolayer systems with 1/r 3 interactions. More- 
over, in the Molecular Dynamics done by S. Lin and 
co-workers 58], the coupling constants for both transi- 
tions, Fluid/Hexatic and Hexatic/Solid, that they found 
in their analysis agree very well with the ones found in the 
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(b) 



FIG. 2: (Color online) Snapshots of monolayer with Voronoi 
and Delaunay constructions near the Fluid/Hexatic transi- 
tion, (a) The parameters are n = 0.71, p = 0.3 and T = 191.7 
; for this thermodynamical state, the MC average of the Bond 
orientational order parameters is < fa >= 0.67 ± 0.03. (b) 
n = 7.30, p = 0.3, T = 157.9 and < fa >= 0.61 ±0.07. White 
cells have 6 sides ; green : 5 sides ; red : 7 sides ; yellow : 4 
sides and blue : 8 or more sides. 



present work done with Monte Carlo computations and 
Ewald sums. Therefore, it seems reasonable to consider 
that the real interaction between the superparamagnetic 
colloids in an external magnetic field deviate slightly form 
a true 1 /r 3 potential ; one should note also that the ex- 
periments and the simulations are in agreement with the 
small amplitude Ar ~ 3 in which the Hexatic phase can 
be observed for n = 3. 

With the empirical curves given by Eqs. (l9l8p and the nu- 
merical value of the Tj-parameters, we found p[ F ^ H ^ > 
p ( c H/S) for 4.1 < n < 127 with T = 1.0 and Q = 14. This 
indicates that the hexatic phase would disappear in this 
range of n and that the KTHNY transitions are replaced 
by a first order phase transition (25j . 
On FIG12 we show snapshots of monolayer systems in 
Hexatic and Hexagonal phases near melting (Runs B). 
For n = 0.71 (FIG^a)), the disclinations are bound in 
pairs (equivalent to isolated dislocations) and quartets 
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FIG. 3: (Color online) Dependence of the structure with tem- 
perature for n = 0.71 and p = 0.3. (a) Center to center 
correlation functions. The fit of the tail is represented for 
T — 1.01. (b) Bond orientational correlation functions. The 
fit is given by Eq.(|3]) and it is represented for T = 1.01; inset 
in (a) shows the agreement between the fit and MC results 
for g&(s). (c) Bond orientational order parameter as function 
of the temperature, < fa > MC is obtained with MC compu- 
tations and < fa > g is the values obtained by fitting the tails 
of ge (s). Inset : Histograms of fa for several temperatures. 



and they are more or less uniformly distributed in the 
monolayer, this is consistent with the KTHNY theory 
[26l [28| and the system is in the fluid phase very near 
to the hexatic transition for the thermodynamic param- 
eters given in the caption of FIGj2ja). For n — 7.30 
(FIGj2Jb)), the defects form small lines and aggregate 
; these arrangements of the defects may be considered 
as seeds for grain boundaries in a first order melting 
[TTI . [25| . Qualitatively, the range of the interactions can 
help us to understand the differences found in the distri- 
bution of defects. In one hand, for long ranged interac- 
tions (n = 0.71), each particle interacts with all particles, 
thus it is energetically unfavourable for the defects to ag- 
gregate in lines. On the other hand, for much shorter 
interactions (n — 7.30), heterogeneous distributions of 
defects are not averaged out by the range of interactions 
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FIG. 4: (Color online) Dependences on temperature and Kosterlitz-Thouless singuralities for n — 0.71 and p = 0.3. (a) Finite 
size analysis of BOOP ; t = (T - T F/H (N))/T F/H (N), the black straight lines correspond to In < (j> e >~ -0.412 + 0.37 | t 1/2 
for t < 0. (b) Finite size analysis of the susceptibility. Inset : T F / H as function of the number of particles, (c) Correlation 
length £ (T) as function of temperature in the fluid phase. Inset : exponent rj(T) in the hexatic phase 0.9 < T < 1.0 and in 
solid phase T < 0.9. (d) Bond orientational correlation lengths £e(T) and £+(T) as function of temperature in the fluid phase. 
Inset : stiffness Ka(T) in the hexatic phase computed as Ka(T) = 18kBT/nrja(T). 



; on contrary, the quite short range of the interaction 
favours the growth of grain boundaries. 



IV. MONTE CARLO COMPUTATIONS AND 
ANALYSIS. 

In FIG [TJ(b), we represent BOOP for Runs A and B, 
and on FIGEJc) they are represented as function of tem- 
perature for n = 0.71 and p = 0.3. 

The insets (c,d) in FIG[T](b) show straight line fits of 
the fourth order cumulants in the region of transition for 
several size of subsystems. For Runs A, FIG[TJ(c), the 
straight line does not cross in a single point ; this argues 
in favour of a first order transition for the melting of the 
hexagonal 2D solid phase [Hj]. For Runs B, FIGEQ(d), 
the shift in the crossing of the straight line fits of the cu- 



mulants is less apparent ; in consideration with the other 
results described in this paper (and done by others in 
previous studies), it seems plausible that the transition 
is weakly first order or located quite close to a kind of 
multicritical point separating the KTHNY theory from a 
first order phase transition 30]. 

On Figs|3][a,b), the red curves are fits obtained with 
Eqs. (|3l4[) . The values of < 4>§ > g , as well as the Monte- 
Carlo averages < <j)6 > mc are reported on FigJ3Jc). 
Both evaluations of the bond orientational order param- 
eter agree extremely well in hexatic and solid hexagonal 
phases, but a significant difference is seen in the fluid 
phases. This difference in the fluid phase is well under- 
stood since bond orientational order in the fluid phase 
is short ranged and < ^ > uc is the average over all lo- 
cal order parameters, while < cf>Q > g is the asymptotic 
value of ge(s). In the solid phase T < 0.90, the cen- 
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FIG. 5: (Color online) Histograms of energies and bond ori- 
entational order parameter <j>6 for several temperatures with 
n — 1.65 and p = 0.109. The choice of n corresponds to uh 
in FigfT] (a-inset) . For each histogram shown, the value of the 
temperature is given in legend. The histograms are computed 
from the trajectories of the MC computations with 9 x 10 
and 2 x 10 6 MC-cycles. These histograms of energies (tra- 
jectories) are used in the multiple histogram method to build 
the best estimate of the density of state, (a) Histograms of 
energy per particle, (b) Histograms of the order parameter 

06- 



tcr to center correlation functions g(s) are exactly those 
of a hexagonal lattice with thermal fluctuations and for 
s > 10 a, the bond orientational correlation functions 
ge(s) —< (f>6 >g=< 06 >mc within the statistical fluc- 
tuations, this indicates a long range of the bond orien- 
tational order. The dependence on temperature for the 
histograms of </>6 (Inset in FIG[3lc) - n = 0.71) is similar 
to the dependence found experimentally in superparam- 
agnetic colloids (n = 3) confined at an air/ water interface 
(cf. FIG.6 of ref.[H|]). For T = 1.00 and T = 1.01, the 
histograms of <f)§ are very broad, it corresponds to the 
KT-singularity of \6 ] similar behaviors are observed on 
histograms computed for n — 1.65 given in FIGEJb). 
On FIG|H we summarize the analysis on the done on 
Kosterlitz-Thouless essential singularities to obtain the 
transition temperatures T f /h and Tjj/g for a set of 
Monte Carlo computations (n = 0.71 and p — 0.3). All 
singularities showed on FIG 0] are plotted with v = 1/2 



; however, the numerical results given in FIG|4]may also 
be quite well represented with 0.35 < v < 0.55. In 
KTHNY theory, the essential singularities as in Eq.([7]) 
occur with v =1/2 (F/H) or with v = 0.369 634... (H/S) 
H M, HI H! ill (or with i/ = 2/5 The values of 

the transition temperatures obtained by using v — 1/2 or 
v are in the uncertainties ; for instance, for X6 in FIG|4jb) 
with v = 1/2 one finds T F / H = 0.9783 and with v one 
has Tp/H — 0.9845, and both fits are equally good. 
At the Fluid/Hexatic transition, the KTHNY theory 
states that v = 1/2, this is fully consistent with the 
Monte Carlo results obtained in the present work ; how- 
ever, we are not able yet to estimate with a convenient 
accuracy the value of the exponent v in Eq.([7]) at the 
Hexatic/Solid transition and verify, as stated in KTHNY 
theory, that v{H/S) = v ~ 0.369 634.... 
On FIG[4ja,b), we give finite size analysis for BOOP and 
susceptibility represented as function of the scaling vari- 
able t=(T- T F/H {N))/Tp /H (N) and, on FIGgfod), 
we represent as function of the temperature £(T), £+(T), 
£ 6 (T), n(T) and rj 6 (T) obtained by fitting the tail of the 
correlation functions. The numerical values of T F / H (N) 
used in the definition of the scaling variable t are ob- 
tained from Eq. ([TJ) for each subsystem size. 
The collapse of data for different subsystem sizes on 
the universal Kosterlitz-Thouless essential singularity for 
both the bond orientational order parameters and suscep- 
tibilities, FIG|4]Ja,b), and also the dependence on tem- 
perature of correlation lengths, FIG0|c,d), are strong 
arguments in favour of the KTHNY theory for n — 0.71 
and p = 0.3. Same results are obtained for Runs A and 
B (cf. TABLE QJ for n < 2 and also for computations 
with n = 1.48, 1.65 and 3. 

For n — 0.71, p — 0.3 and T = 196, compiling the results 
from finite-size analysis and KT-singularities, we found 
that the Fluid/Hcxatic transition is located at temper- 
ature T F/H = 0.986 ± 0.008. From the KT-singularities 
found for the stiffness K A (T) (inset FIGHJd)) and £+(T) 
(blue line on FIGHfd)), we found the temperature of the 
Hexatic/Solid transition at T H/S = 0.913 ± 0.011. From 
the singularity found for Ka{T), we compute t]q(T h / s ) — 
0.385, that is different from the value 1/4 predicted in the 
KTHNY theory ; we may interpret this disagreement as 
it is done for Coulomb interaction [loj . by arguing that 
the critical exponents depend on the range of the inter- 
action [71], [lH . More precisely, we consider that the peri- 
odic boundary conditions, which introduce an arbitrary 
long range positional order between particles, and, more 
fundamentally, the dependence of critical ex pon ents with 
the range of the interaction potentials [lol l7ll. |72T | may 
explain the difference between the KTHNY predictions 
in the Hexatic phase and the observed critical behavior 
of £(T) and values of r)(T) and r] 6 (T) reported on insets 
of FIG0|c,d). A complete study of the dependence of 
the critical exponents with n is needed to understand 
precisely the critical behavior of the 71-OCP monolayer 
systems and modifications induced by the range of inter- 
actions in the KTHNY predictions ; such a study, would 
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FIG. 6: (Color online) Interpolations of thermodynamic observables with the multiple histogram method as function of the 
temperature for n = 1.65 and p = 0.109. All lines, but the KT singularity in (d), are computed with the multiple histogram 
method and symbol are averages computed with Monte Carlo simulations. Data in black are for N = 4096 and curves in blue 
to N = 2025 ; for clarity, the Monte Carlo averages for N = 2025 are not represented on the figures, (a) Average energy 
per particle as function of T (Inset : specific heat as function of T). (b) Bond orientational order parameter <j>$ and 4>\2- (c) 
Interpolations of the fractions of m-coordinated particles computed with Voronoi constructions, x m are defined as < N m > /N 
where N m is the number of particles with exactly m neighbors for a given configuration. The quantity 1 — x§ is the fraction 
of defects, (d) Susceptibility of the order parameter (f>e- The KT singularity is obtained by fitting the multiple histogram 
interpolation of the susceptibility by using Eq.® with v = 1/2, the fitted values are Xq — 0.322, a ~ 0.823 and T F / H ~ 0.868. 



be of great interest in the understanding of the melting 
in two dimensions. 

As shown in the following, the Kosterlitz-Thouless essen- 
tial singularity can be very well represented with the help 
of the multiple histogram method (MHM) (6f| [66| . On 
FIGfSl we represent the energy and the order parameter 
histograms for the temperatures used in the multiple his- 
togram reweighting method the system with N = 4096 
particles for n — 1.65 and p — 0.109. The histograms 
shown on FIGf5] are computed from the trajectories of 
the Monte Carlo simulations. The overall qualitative fea- 
tures of the order parameter histograms for both n — 1.65 
(FIGGJb)) and n = 0.71 (Inset of FIG^c)) are similar 
and in agreement with the ones found for superparamag- 
netic colloids confined at an air/water interface |l3[ (see 



also FIG. 2 of ref.[18(). It is also worthwhile to note that, 
close to the KT-singularity (T = 0.90-0.88) the order pa- 
rameter histograms broaden strongly and tend to become 
quite flat between for <f>Q between 0.1 and 0.6, making the 
susceptibility x& violently increasing. From the MC tra- 
jectories used to compute the histograms of FIGfSl we 
apply the multiple histogram reweighting technique to 
obtain interpolations of the thermodynamic observables 
at temperatures other than the ones used in the Monte 
Carlo computations. These interpolations are shown on 
FIGJ6] for systems with 4096 particles (black curves and 
symbols) and 2025 particles in blue. 
On FIG HJa), we represent the average energy per particle 
as function of T ; the variation of the energy above Tp/H 
is due to the gradual dissociation of disclination pairs 
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FIG. 7: (Color online) Center to center correlation functions 
g pq (s) for pair of particles with p neighbors and q neighbors for 
n = 1.65 and p — 0.109. For the <?57(s) correlation function 
(top), at T = 0.83, the first two peaks are approximately 
located at s — 3.1 a (height : ~ 110) and 6.1 a < s < 6.6 a 
(height : ~ 7). The first two peaks in the gr? (s) correlation 
function (bottom) are located at s — 4.2 a (height : ~ 45) 
and 6.0 a < s < 6.6 a (height : ~ 12). 



(see also FIG|SIc)). In the inset we represent the specific 
heat computed with the Multiple Histogram Method for 
both system size N = 4096 and N = 2025, the transition 
temperatures Tp/jj and Tjj/s found in KT-singularities 
are represented as vertical lines. The peak in the spe- 
cific heat near the Fluid/Hexatic transition occurs above 
T F / H and it is also due to gradual dissociation of disclina- 
tion pairs in agreement with the KTNHY theory [24l [H[ 
; this peak is similar to the one found because of the 
vortex contribution in the XY-model according to which 
the height, position and shape of the maximun are not 
universal (2|| H3|. On the system with 4096 particles, a 
small bump is seen also on the specific heat just before 
Tjj/g, it would correspond to the dissociation of disloca- 
tion pairs. 

On FIGEfb), we represent the bond orientational order 
parameters (j> 6 (N = 4096 and 2025) and <f> 12 (N = 4096) 
and on FIGlSJd), the susceptibility \6 is shown for both 
system sizes. The effects of the finite size of the sys- 
tem (and periodic boundary conditions) can be seen on 



both figures and these finite size effects are consistent 
with the theory j&ij ]. The maximum found for the sus- 
ceptibility with the multiple histogram method is very 
slightly shifted to the higher temperature with decreas- 
ing the system size. More precisely, for the system with 
N = 4096 particles the maximum is located at T = 0.886 
and xe( T = 0.886 ;4096) = 114.3 while for the system 
with N = 2025 particles the maximum is at T = 0.888 
and xe(r = 0.888 ;2025) = 54.6, this agrees well with 
the finite size analysis done on subsystems (cf. FIG|U[b) 
- symbols for N = 4096 and < N >= 2048). 
On FigJSJd), the red curve is the KT-singularity, Eq.©, 
obtained by fitting the multiple histogram interpolation 
of the susceptibility, the agreement is very good both 
with the interpolation done with the multiple histogram 
reweighting and with the finite size analysis on subsys- 
tems. 

On FIGHJc), the fraction of m-coordinated particles with 
to = 5 and 7 and the fraction of defects, defined as 1 — xq, 
are shown. On all the temperature range used in the mul- 
tiple histogram interpolation, one has x$ = xr with an 
excellent accuracy in the solid and hexatic phase and we 
found 1 — xg =2:5+2:7. In the fluid, there are few parti- 
cles with to = 4, 8 and 9 neighbors, however, the equality 
X5 = xj still holds to about 1% accuracy. The fractions 
of TO-coordinated particles with to = 4, 8 and 9 are in 
the ordered phases (hexatic and solid) : X4 < 4 x 10~ 5 
; ig < 2 x 10~ 5 and xg < 10 -9 ; and, in the fluid phase 
: x 4 ~ 7 x 10~ 4 , x s ^ 10~ 3 and 2:9 < 10~ 6 . From the 
fraction of defects in the Solid phase, we may estimate 
the core energy E c at low temperature by using 

1 - s 6 = X d °exp f- 2 25 c (n)^) (11) 

to model the defect fraction as an Arrhenius law in the 
low temperature limit [H, 0, 0, EH . 
With Eq.pj, we found : for n = 1.48, E c = 4.38 ; 
n = 1.65, E c = 5.02 and for n = 3.00, E c = 3.46. The 
change from the KTHNY melting to the grain boundary 
induced melting is predicted to occur when E c < 2.84 
[H [H [13, E(f ; therefore, the results found for E c 
from the analysis of the fraction of defects support the 
KTHNY melting according to the Chui-Saito criterium 
[25|, HH (see also next section). 

The center-to-center correlation functions 557(5) and 
577 (s) are represented on FIG(7]for the Fluid and Hexatic 
phases and near the Hexatic/Solid transition (T — 0.83). 
As explained in Sect. II, in Solid and Hexatic phases the 
fraction of defects is small, thus large statistical fluctua- 
tions are found in (757 (s) and 577(5) for s > 10<t ; for in- 
stance, at T = 0.83, one has x 5 = x 7 = (1.3±0.3) x 10~ 2 . 
The first peak in the 557(5) gives an estimate of the 
strength of the bound between the 5 and 7-coordinatcd 
particles (positive and negative disclinations) ; at T = 
0.88 (below T F / H - cf. the KT-singularity of x&)i the 
first peak is located at s = 3.1 a and the height is about 
52 and, just above T F / H , at T = 0.89, the height of the 
first peak falls to 14.2. In the Fluid phase (T = 1.00), 
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FIG. 8: (Color online) Topological defects in the Hexatic (near Solid) phase ; n = 1.65, T = 0.83 and p = 0.109. (a,b,c) Quite 
frequent defects ; (d,e,f) rare defects. In all snapshots, the diameter of the coloured disks is a(= 1). The Delaunay construction 
is represented only for particles that do not have 6 neighbors. The Voronoi cells are represented for all particles. The color 
code for the particles is defined by the number of neighbors : 4 neighbors (yellow) , 5 (green) , 7 (red) , 8 (blue) . The small black 
disks inside white Voronoi cells are the point particles with 6 neighbors. 




FIG. 9: (Color online) Topological defects in the Hexatic (near Fluid) phase ; n = 1.65, T = 0.88 and p — 0.109. 



the first peak is located at s — 3.0 a and its height is 
about 8.5. 

Similarly, the height of the first peak in the 377(5) gives 
an estimate of the strength of the bound between dis- 
locations (see FIGH (a) and (d)). At T = 0.83, in the 
Hexatic phase very close to the Hexatic-Solid transition, 
the first peak is located at s — 4.2 a and its height is 
about 45.1, it falls to 17.5, in the Hexatic phase close to 
the Fluid-Hexatic transition (T = 0.88) and to 3.7 in the 
Fluid phase at T = 0.89. 

For both correlation function 357(a) and 377(3) (and also 



<?55(s) - not shown), the locations of the first peak are 
extremely well represented with the topological defect 
shown on FIG^a), that corresponds to a pair of dis- 
locations or a quartet of disclinations. Thus, the anal- 
ysis of the behaviour of the first peaks of the correla- 
tions functions 357(3) and 377(a) with an increase of the 
temperature supports the view of the KTHNY theory 
that describes the melting from Solid to Hexatic as a 
dislocation-unbinding transition, followed by the Hex- 
atic to Fluid transition stemming from the disclination- 
unbinding transition (24[. 
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More generally, the first three peaks of 357(5) and 377(a) 
are very well located by taking into account only the de- 
fects of FIGE] (a-c). The defects shown in FIG® (b) 
and (c) are quite frequent in the Hexatic phase, how- 
ever, as it can be seen on the snapshot given in FIG[2ja), 
other topological arrangements of disclinations and dis- 
locations with irregular shapes can also be found, some 
are given in FIGJH1 The defects shown on FIG|5](d-f) are 
very rare. 



V. DISCUSSION AND PERSPECTIVES 

In the KTHNY theory, the melting of the hexago- 
nal solid in two dimensions is described by two continu- 
ous transitions. From the hexagonal Solid phase to the 
isotropic Liquid phase, the first transition is induced by 
the unbinding of the dislocations and it corresponds to 
the Solid-Hexatic transition, the second transition is the 
disclination-unbinding transition from the critical Hex- 
atic phase to the isotropic liquid [24], HU, H3| ■ 
From the continuum elasticity theory, the dislocation- 
unbinding transition is described with the Hamiltonian 

m 



come large [2J| (see also FIG|2](a) an d (b)). Therefore, 
the core energy E c can be used to define a criterium for 
the dislocation-unbinding transition : for E c < 2.84, the 
transition is (weakly) first order and for E c > 2.84 it fol- 
lows KTHNY theory. 

In experiments and in computer simulations, we may ob- 
tain an estimate of E c from the fraction of defects in the 
solid phase [12|,|40|,|43|,|44j by using an Arrhenius law as in 
Eq. flTTj) . The estimates of E c done in the previous section 
for the n-OCP monolayer show first that, for n < 3, the 
dislocation-unbinding transition follows KTHNY theory 
according to the Chui-Saito criterium, and, second, that 
E c depends on n. In the present work, the estimates of 
E c are given for only three values of the power n ; this 
does not permit yet to obtain the dependence of E c on n 
and to improve the Chui-Saito criterium. 
A dislocation is equivalent to a pair of disclinations 
strongly bounded. As shown by B.I. Halperin and D.R. 
Nelson [28j . the disclination-unbinding transition can be 
described by using an isomorphism with the XY model 
[3lT ] and the Coulomb gas in two dimensions 29], with 
the disclinations represented as vortices or charges. The 
Hamiltonian for the disclination-unbinding transition is 



K 

8^ 



E 



b(r) ■ (r-r')b(r') ■ (r - r') 



r — r 



1 12 



6(r) • b(r') In |r T 1 (12) 



+£ c £|6(r) 



with b(r) the Burgers vector of a dislocation located at 
r, K the Young's modulus, a c the core radius that de- 
fines the distance outside of which the dislocation energy 
obey the logarithmic form of the previous equation and 
E c is the core energy that gives the contribution to en- 
ergy inside the core radius. 

The core energy E c has a crucial influence on the mech- 
anisms of melting in two dimensions [25l |32| . 1 3 31 ] . Renor- 
malization group analysis of the KTHNY theory uses 
the fugacity of defects, y — exp(—E c ), in the recursive 
relations [26], [4(j. In the analytical description of the 
grain-boundary theory of melting 25], S.T. Chui predicts 
that the melting is weakly first order when E c < 2.84 
and becomes strongly first order as E c > 2.84. How- 
ever, Monte- Carlo simulations using the Hamiltonian 
of Eq. (IT2l done by Y. Saito [H, [33| showed that the 
dislocation-unbinding transition is well described by the 
KTHNY theory for large E c while it became first order 
for small value of E c , because of the nucleation of grain 
boundary loops. In the grain-boundary induced melting 
theory, it is assumed that the distance between dislo- 
cations making a grain boundary is small between grain 
boundary, this approximation breaks down for large value 
of E c for which the distance between dislocation pairs be- 



36 



s(r)s(r') In 



+£cd^s(r) 2 



O-cd 



(13) 



with Ka the Frank constant (the coupling constant re- 
lated to the distortion of the bond-angle field), s(r) the 
charge of the disclination located at r, a c d is the disclina- 
tion core size and E c d the disclination core energy. In the 
isomorphism with the Coulomb gas, the charges assigned 
to disclinations can be chosen as s(r) = +1 for a particle 
with seven nearest neighbors (Voronoi cells and parti- 
cles in red in FIGSI2J [3] and and —1 for five nearest 
neighbors (represented in green). The analysis of the first 
peaks of the correlation functions <?57(s) and <777(s) (cf. 
FIGfT]) done in the previous section is in excellent qualita- 
tive agreement with the disclination-unbinding transition 
for the Hexatic and Fluid phases. 

The results of the Kosterlitz-Thouless theory for the 
XY model and Coulomb gas in two dimensions can be 
used for the disclination-unbinding transition, the Hex- 
atic/Fluid transition. In the XY model, the correlation 
length has an essential singularity as in Eq.([7]) |3l|. All 
quantities related to the correlation length exhibit the 
same KT essential singularity at the transition, for the 
ri-OCP monolayer, KT-singularities are found for : the 
bond orientational order parameters 4>§ (FIG0](a)), the 
susceptibility \6 (FIGHJb) and FIGJBfd)) and correla- 
tions lengths f (T) and £ 6 (T) (FIGgfod)). Tn e transition 
temperatures T F / H or coupling constants T F / H found by 
fitting the KT-singularities of these quantities agree well 
and permit to locate the Hexatic/Fluid transition for sev- 
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eral values of the power n of inverse power law interac- 
tion of Eq.([T]). These results are reported on TABLE 
HI As shown also in the previous section, the Kosterlitz- 
Thouless essential singularities can be very well studied 
with the multiple histogram method ; additional compu- 
tations with different sizes and others values of n are on- 
going. These results for n < 3 are in excellent agreement 
with the KTNHY theory for the disclination-unbinding 
transition, they agree also with previous experimental 
[I SOS Ql and numerical HHSH studies. 
In the present work, as n is a continuous parameter in 
the IPL potential, we have been able to obtain empir- 
ical transition curves for both Solid/Hexatic and Hex- 
atic/Fluid transition for n < 3. These curves are given by 
Eqs. (|8l9[) in Sect. Ill and they are represented on FIGfjTa) 
in the plane (n,p). This phase diagram is extremely 
useful to locate the Hexatic phase for a given n and to 
choose the range of numerical values for the surface den- 
sity p, the temperature T or the coupling constant V in 
numerical computations that one has to use to observe 
the Hexatic phase, this is necessary to study in detail the 
properties of the Hexatic phase and of the Hexatic/ Solid 
transition since the range in the coupling constant V for 
which one can observe an Hexatic phase is very thin for 
almost all values of n < 3. It is the reason for which, 
we consider the results obtained in this work as prelim- 
inary results for a longer study of the Hexatic phase in 
the n-OCP monolayers. Obviously, the accuracy of the 
transition curves reported in Sect. Ill, and the I\ param- 
eters given in TABLE HI will be improved by adding data 
obtained with other values of n. 

For n > 3, we found that the transition curves cross at 
n ~ 4, making pc > Pc \ and computations done at 
constant surface density and temperature (Runs A and 
B - FIGQJ indicate that the melting would not follow 
the KTHNY theory for n > 4. The transition curves 
given in Sect. Ill do not permit to reproduce to a cor- 
rect accuracy the data obtained in Runs A and B for 
n ~ 12.7 and n ~ 7.5. However, the finite size analysis of 
the fourth order cumulants using subsystems and some 
double peaked histograms of the bond orientational order 
parameters observed for n = 12.75 — 12.95 may support 
the assumption that the melting becomes a first order 
transition. Nevertheless, we are not yet able to reach a 
definitive conclusion about the nature of the melting for 



n > 4 and more computations are needed for these values 
of n. 

The recent results on the hard disks systems described 
by E.P. Bernard and W. Krauth [38|], show that for 
n — > oo the melting would be described by a first order 
Fluid/Hexatic transition followed by a second order Hex- 
atic/Solid transition. These results do not agree neither 
with the KTHNY theory, nor with the grain-boundary 
induced melting proposed by Chui. However, these re- 
sults might find a nice interpretation with the help of 
the theory of the melting based on isostructural critical 
points [34| - l36j in which a first order transition between 
the Fluid and Hexatic phases is plausible [36] ; this the- 
ory has also to be taken into account for studies with 
large but finite n. 

The empirical transition curves for n < 3 and the excel- 
lent qualitative and semi-quantitative agreements with 
the KTHNY theory found for the melting of n-OCP 
monolayers with n < 3 are the main results of the present 
work. 



Acknowledgments 

I am very grateful to Dominique Levesque, Jean- 
Michel Caillol, Jean- Jacques Weis, Emmanuel Trizac, 
Gerhard Kahl and Moritz Antlanger for many inter- 
esting discussions. The author acknowledges also the 
computation facilities (iDataPlex - IBM) provided by 
Direction Informatique of Universite Paris-Sud. 



Appendix A: Ewald sums for Inverse Power Law 
interactions. 

In this appendix, we give the Ewald formulas used to 
compute the energy of the monolayer. A complete deriva- 
tion of the Ewald method starting from the interaction 
potential of Eq. (TT]) and for any continuous values of n can 
be found in refelgllc}. 

For n ^ 2, the energy of a configuration of the n-OCP 
monolayer is given by 
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with Si the location of the point particle in the mono- the system, G are the wave vectors of the reciprocal lat- 
layer, S n a symbolic notation for the periodic images of 
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tice associated to the periodic lattice S n and a is the 
Ewald parameter [49| . For n = 2, the energy is given by 
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with E±(x) = —Ei(—x) the exponential integral. 
In Eq. flAip . T(a;x) is the complementary incomplete 
gamma function, when a < (i.e. n > 2) the function 
has to be computed by using the recursive relation 



T(a + 1; x) = ar(a; x) + x a e~ 



(A3) 
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